SoSe2021
Folienübersicht
Lineare Modelle in R
Anwendungsbeispiel iris
lm()
Wir modelliere die Kronblattlänge als Funktion der Kelchblattlänge:
mod <- lm(formula = Petal.Length ~ Sepal.Length, data = iris) mod
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)
Coefficients:
(Intercept) Sepal.Length
-7.101 1.858
coef(mod)
(Intercept) Sepal.Length -7.101443 1.858433
Wie lautet die vollständige lineare Regressionsgleichung?
Populations- vs. Stichprobengleichung
1
\[y = a + b*x\]
- Die häufigste Art, die Gleichung zu schreiben ist die Stichprobengleichung für die erklärenden y-Werte!
Zur Erinnerung
- Parameter = eine Messgröße, die eine POPULATION beschreibt (z.B. \(\mu, \sigma^{2}\))
- Statistik = eine Messgröße, die eine STICHPROBE (‘sample’) beschreibt (z.B. \(\bar{x}, s^2\))
Populations- vs. Stichprobengleichung
2
Populationsgleichung:
- für beobachtete y-Werte: \(y_{i} = \alpha + \beta*x_{i} + \epsilon_{i}\)
- für vorhergesagte y-Werte: \(\mu(y_{i}) = \alpha + \beta*x_{i}\)
Stichprobengleichung:
- für beobachtete y-Werte: \(y_{i} = a + b*x_{i} + \epsilon_{i}\)
- für vorhergesagte y-Werte: \(\hat{y}_i = a + b*x_{i}\)
Die Stichprobenkennwerte \(a\) und \(b\) repräsentieren die Punktschätzer der Regressionsparameter \(\alpha\) und \(\beta\)!
Anwendungsbeispiel iris
coef(mod)
(Intercept) Sepal.Length -7.101443 1.858433
Für die beobachteten Stichprobenwerte gilt:
\[Petal.Length_i = -7.1 + 1.85*Sepal.Length_i + \epsilon_{i},~~\text{wobei}~~\epsilon_i ~N(\mu, \sigma^2)\]
R rechnet noch viel mehr aus: die summary() Funktion
summary(mod)
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)
Residuals:
Min 1Q Median 3Q Max
-2.47747 -0.59072 -0.00668 0.60484 2.49512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.10144 0.50666 -14.02 <2e-16 ***
Sepal.Length 1.85843 0.08586 21.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared: 0.76, Adjusted R-squared: 0.7583
F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16Wo ist a? Wo ist b? Und was bedeuten die anderen Zahlen…?
Die Genauigkeit der Schätzung
Von der Stichprobe auf die Grundgesamtheit verallgemeinern
- Wie wären die geschätzten Koeffizienten, wenn wir das Experiment wiederholen würden?
- Wo liegen die geschätzten Koeffizienten für die gesamte Population?
- Hier braucht es die Berechnung der Standardfehler und Hypothesentests
- Zur Erinnerung: Die Genauigkeit der Punktschätzung sollte immer durch den Standardfehler des Schätzers angegeben werden → je kleiner der Standardschätzfehler, desto größer ist die Effizienz des Schätzers.
Reproduzierbarkeit geschätzter Koeffizienten
Hohe Varianz
Stelle Dir vor, jeder Datenpunkt repräsentiert einen individuellen Fisch und alle Datenpunkte zusammen repräsentieren eine Population.
Wenn Du jetzt mehrere Zufallsproben von 10 Fischen nimmst, wie ähnlich wären sich dann die geschätzten Koeffizienten?
Reproduzierbarkeit bei starker Streuung
N = 10
Rate: Welche sind sich viel ähnlicher und welche sind näher an den wahren Populationskoeffizienten dran?
Reproduzierbarkeit bei starker Streuung
N = 10
‘Wahre’ Parameter: \(\alpha\) = 11.27 und \(\beta\) = 0.41
Reproduzierbarkeit bei starker Streuung
N = 100
Lass uns jetzt die Probengröße auf 100 Fische steigern!
Reproduzierbarkeit bei starker Streuung
N = 100
‘Wahre’ Parameter: \(\alpha\) = 11.27 und \(\beta\) = 0.41
Reproduzierbarkeit bei schwacher Streuung
N = 10
Wie ähnlich sind sich die geschätzten Koeffizienten bei schwacher Streuung und einem Stichprobenumfang von 10 Fischen?
Reproduzierbarkeit bei schwacher Streuung
N = 10
‘Wahre’ Parameter: \(\alpha\) = 10.04 und \(\beta\) = 0.5
Reproduzierbarkeit bei schwacher Streuung
N = 100
Lass uns jetzt die Probengröße auf 100 Fische steigern!
Reproduzierbarkeit bei schwacher Streuung
N = 100
‘Wahre’ Parameter: \(\alpha\) = 10.04 und \(\beta\) = 0.5
Maß der Genauigkeit der geschätzten Koeffizienten: \(s_\bar{a}\) und \(s_\bar{b}\)
Die Unbestimmtheit der geschätzten Steigung und Achsenabschnitt
- steigt mit zunehmender Varianz und abnehmender Probengröße
- ist größer wenn die Spanne der Werte klein ist (gemessen als \(SS_{X}\))
Standardfehler des Achsenabschnitts und der Steigung:
\[s_\bar{a}=\sqrt{\frac{MS_{Residual}\sum x_i^2}{n*SS_x}}=\sqrt{\frac{\frac{SS_{Residual}}{df}\sum x_i^2}{n*SS_x}}=\sqrt{\frac{\frac{\sum(y_i-\hat{y_i})}{n-2}\sum x_i^2}{n*\sum(x_i-\bar{x})^2)}}\]
\[s_\bar{b}=\sqrt{\frac{MS_{Residual}}{SS_x}}=\sqrt{\frac{\frac{SS_{Residual}}{df}}{SS_x}}=\sqrt{\frac{\frac{\sum(y_i-\hat{y_i})}{n-2}}{\sum(x_i-\bar{x})^2}}\]
Zur Erinnerung: Varianzzerlegung in linearen Modellen
Die Gesamtstreuung (\(SS_Y\) oder \(SS_{Gesamt}\)) kann in 2 Komponenten zerlegt werden:
erklärbare Streuung = \(SS_{Regression}\) + nicht erklärbare Streuung = \(SS_{Residuen}\)
Was sind die Standardfehler im Kronblattmodell?
summary(mod)
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)
Residuals:
Min 1Q Median 3Q Max
-2.47747 -0.59072 -0.00668 0.60484 2.49512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.10144 0.50666 -14.02 <2e-16 ***
Sepal.Length 1.85843 0.08586 21.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared: 0.76, Adjusted R-squared: 0.7583
F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16→ Punktschätzung für den Achsenabschnitt: -7.10144 ± 0.50666
→ Punktschätzung für die Steigung: 1.85843 ± 0.08586$
Konfidenzintervalle der ‘wahren’ Koeffizienten
confint()
- Konfidenzintervalle kann man nicht nur für Mittelwerte oder Varianzen berechnen, auch für Regressionskoeffizienten.
- Berechnung basiert auf der t-Verteilung:
Manuelle Berechnung
a <- -7.10144
b <- 1.85843
ci_a <- 0.50666 * qt(p = 0.975, df=148)
ci_b <- 0.08586 * qt(p = 0.975, df=148)
mat <- matrix(
c(a-ci_a, a+ci_a, b-ci_b, b+ci_b),
nrow = 2, byrow = TRUE)
rownames(mat) <- c("(Intercept)", "Sepal.Length")
colnames(mat) <- c("2.5%", "97.5%")
mat
2.5% 97.5% (Intercept) -8.102662 -6.100218 Sepal.Length 1.688760 2.028100
Mit build-in Funktion
confint(mod, level = 0.95)
2.5 % 97.5 % (Intercept) -8.102670 -6.100217 Sepal.Length 1.688772 2.028094
- Mit 95%iger Wahrscheinlichkeit liegt \(\alpha\) zwischen -8.1 und -6.1 und \(\beta\) zwischen 1.7 und 2.0
- Beide KI schließen nicht die Null mit ein → wir können davon ausgehen, \(\alpha \neq 0\) und \(\beta \neq 0\).
Signifikanztest der Koeffizienten
- Die Schätzungen der Koeffizienten basierend auf der Stichprobe weichen immer von Null ab.
- Ob dies Zufall ist oder nicht kann mit einem t-Test jeweils geprüft werden:
- \(H_0:~\alpha = 0\), d.h. die Regressionsgerade geht durch den Ursprung.
- \(H_0:~\beta = 0\), d.h. der Mittelwert ist für die Schätzung ebenso geeignet wie die Regressionsgerade.
- Die Teststatistik t ist einfach nur der Quotient aus dem Schätzwert und seinem Standardfehler und folgt einer t-Verteilung:
- \(t_a = \frac{a}{s_\bar{a}}~~\text{und}~~t_b=\frac{b}{s_\bar{b}}\)
Teststatistiken
(t_a <- -7.10144 / 0.50666)
[1] -14.01618
(t_b <- 1.85843 / 0.08586)
[1] 21.64489
Wahrscheinlichkeiten
pt(t_a, df = 148)
[1] 3.065741e-29
pt(t_b, df = 148, lower.tail = FALSE)
[1] 5.224047e-48
Testergebnisse im summary() output
summary(mod)
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)
Residuals:
Min 1Q Median 3Q Max
-2.47747 -0.59072 -0.00668 0.60484 2.49512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.10144 0.50666 -14.02 <2e-16 ***
Sepal.Length 1.85843 0.08586 21.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared: 0.76, Adjusted R-squared: 0.7583
F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16
Interpretation der p-Werte
Wenn ein p-Wert sehr niedrig ist ( < 0.05), bedeutet das, dass
- der wahre \(\alpha\) und \(\beta\) Wert nicht 0 ist.
- Du eine ungewöhnliche Stichprobe hast.
- Deine Annahmen für das Modell falsch sind → daher immer anschließend eine Modelldiagnostik durchführen!
Wie gut ist das Modell?
Zur Erinnerung: das Bestimmtheitsmaß \(R^2\)
Die statistische Kenngröße \(R^2\) ist ein Maß für den Anteil der Variabilität in y, welcher durch das lineare Modell erklärt werden kann:
\[R^2 = \frac{SS_{Regression}}{SS_{Gesamt}} = \frac{\sum(\hat{y}_i-\bar{y})^2}{\sum(y_i-\bar{y})^2}\] \[1 \geq R^2 \geq 0\]
- Werte nahe 1 weisen auf eine sehr gute Beschreibung der Daten.
- Werte nahe 0 weisen auf eine sehr hohe Streuung der Daten um die Gerade an
- Bei \(R^2=1\) lägen alle Datenwerte exakt auf der Geraden.
- In der Ökologie wird bereits ein \(R^2\) von 0.5 bzw. 50% als sehr gut erachtet.
- \(R^2=1\) ist auch das Quadrat des Korrelationskoeffizienten \(r\)
Das Bestimmtheitsmaß \(R^2\) im summary() output
summary(mod)
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)
Residuals:
Min 1Q Median 3Q Max
-2.47747 -0.59072 -0.00668 0.60484 2.49512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.10144 0.50666 -14.02 <2e-16 ***
Sepal.Length 1.85843 0.08586 21.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared: 0.76, Adjusted R-squared: 0.7583
F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16
Modellvalidierung
Annahmen bzw. Voraussetzungen für lineare Regressionsmodelle
- Unabhängigkeit (am Wichtigsten!)
- Homogenität / homogene Varianz (am Zweitwichtigsten)
- Normalität / Normalverteilung
- Linearität
Annahmen/Voraussetzungen
Unabhängigkeit
- Freiheitgrade inkorrekt (Anzahl der unabhängigen Komponenten) → dadurch Beeinflussung der p-Werte (meist stat. Fehler 1. Art erhöht)
- Gibt
- Abhängigkeit die man “riechen” kann
- Abhängigkeit durch nicht passende Modelle
- Zeitliche oder räumliche Abhängigkeit
Annahmen/Voraussetzungen
Varianzhomogenität
Es wird angenommen, dass
- die Varianz in Y konstant ist (z.B. die Varianz bleibt unverändert trotz größer/kleiner werdender y-Werte)
- nur eine Varianz berechnet werden muss und nicht für jeden x-Wert.
- die Residuen hinsichtlich jeder X-Variable homogen sind.
p
Annahmen/Voraussetzungen
Replikation
Normalität und Varianzhomogenität bei Replikation offensichtlicher
Angenommen Du hast 12 Messwerte für die Variable x und y und willst den (positiven) Einfluss von x quantifizieren:
- → Überprüfung beider Annahmen anhand einer Exploration der y Variable schwierig
- → Was aber, wenn es pro X-Wert 50 Beobachtungen gäbe?
Normalität und Varianzhomogenität bei Replikation
Vergleich
Die jeweils 50 Beobachtung bei x = 2,3,4… sollten normal verteilt sein und die Varianzen zwischen den jeweils 50 Beobachtungen homogen:
Normalität und Varianzhomogenität bei Replikation
Vergleich
Die jeweils 50 Beobachtung bei x = 2,3,4… sollten normal verteilt sein und die Varianzen zwischen den jeweils 50 Beobachtungen homogen:
Annahmen/Voraussetzungen
Linearität
- Wenn eine gerade Linie nicht gut die Beziehung zwischen 2 Variablen beschreiben kann,
- sollten die Variablen transformiert werden für die Anwendung einer linearen Regression
- oder es kann eine polynomiale Regression als Spezialfall der multiplen linearen Regression angewendet werden.
- Die quadratische Funktion z.B. hat drei Parameter und kann verschiedene Typen von Parabeln produzieren: \(y=a+bx+cx^2\)
Wie können die Annahmen überprüft werden?
VOR dem Modellieren
- Varianzhomogenität:
- Streudiagramm mit X vs. Y
- Normalität:
- Histogramm von Y
- Shapiro-Wilk-Test mit Y
- Linearität
- Streudiagramm mit X vs. Y
NACH dem Modellieren
Führe eine Modellvalidierung durch und erstelle die sog. Diagnostikplots mit den Residuen:
- Unabhängigkeit:
- Plotte Residuen vs. jede X-Variable
- Zeitserien: Plotte die ACF Funktion der Residuen, wende den Durbin-Watson-Test auf die Residuen an
- Räumliche Daten: Erstelle ein Variogramm mit den Residuen
- Varianzhomogenität:
- Streudiagramm mit Residuen vs. vorhergesagte Y-Werte
- Normalität:
- Histogramm oder QQ-Plot der Residuen
Überprüfe Ausreißer/einflussreiche Beobachtungen
Leverage / Cook’s Distanz ⇒ weiterer Schritt der Modellvalidierung
- ‘Leverage’ h: ein Werkzeug, dass Extremwerte der erklärenden Variable identifiziert, welche potentiell auch das Regressionsergebnis beeinflussen könnten.
- Cook’s Distanz D: das Maß für zu einflussreiche Datenpunkte. Identifiziert einzelne Beobachtungen, die alle Regressionsparameter beeinflussen. Es berechnet eine Teststatistik D, die die Veränderung aller Regressionsparameter misst, wenn eine Beobachtung ausgeschlossen wird.
- D > 0.5 kritisch
- D > 1.0 sehr kritisch
→ Es ist einfacher den Ausschluss beeinflussender Datenunkte zu rechtfertigen, wenn eine große Cook’s Distanz und eine große Leverage für den Datenpunkt vorliegt.
Überprüfe Ausreißer/einflussreiche Beobachtungen
Leverage / Cook’s Distanz ⇒ weiterer Schritt der Modellvalidierung
p
Standardgrafiken zur Modellvalidierung
par(mfrow = c(2,2)) plot(mod)
Gewöhnliche vs standardisierte Residuen
- Gewöhnliche Residuen (\(e_i\)): beobachtete Werte - vorhergesagte Werte (‘observed – fitted’)
- Standardisierte Residuen (\(e_{stand.}\)): Residuum dividiert durch seine Standardabweichung: \[e_{stand.} = \frac{e_{i}}{\sqrt{MS_{Residual}*(1-h_{i})}}\]
- wobei \(h_{i}\) = leverage der einzelnen Beobachtungen i; \(MS_{Residual}\) = Mittlere Summenquadrate der Residuen
- Standardisierte Residuen werden als normalverteilt angenommen mit einem Mittelwert von 0 und einer Varianz von 1; N(0,1).
- Folglich: große Residuenwerte (>2) bedeuten eine schlechte Anpassung des Regressionsmodels.
Berechnen der Residuen in R
Du kannst beide Arten von Residuen berechnen mit
residuals(model)(oderresid()) vom ‘stats’ Paket → ergibt einen Vektor mit den gewöhnlichen Residuenrstandard(model)vom ‘stats’ Paket → ergibt einen Vektor mit den standardisierten Residuenadd_residuals(data, model, var = "resid")vom ‘modelr’ Paket (in ‘tidyverse’) → fügt die Variable ‘resid’ mit den gewöhnlichen Residuen zu deinem ‘data frame’ hinzu; hilfreich beim ‘piping’ von Funktionsaufrufen!
Berechnen der vorhergesagten Werte in R
Du kannst Dir die vorhergesagten Werte für dein Modell über 3 Funktionen ausgeben lassen:
fit(model)→ gibt einen Vektor mit den angepassten Werte zurück, die vom Modell für die Werte der erklärenden Variablen vorhergesagt werdenpredict(model)von dem internen ‘stats’ Paket → verwendet Informationen aus dem angepassten Modell, um eine Glättungsfunktion für die Darstellung einer Linie im Streudiagramm zu erstellen. Der Funktion kann ein data frame mit Sequenzen von Werten aller X Vaiablen, die im Modell sind, übergeben werden (Argumentnewdata).add_predictions(data, model, var = "pred")von dem ‘modelr’ Paket (in ‘tidyverse’) → fügt die Variable ‘pred’ mit den vorhergesagten Werten zu deinem Datensatz hinzu; hilfreich bei ‘piping’ Operationen!
Weitere Diagramme zur Modellvalidierung
Modellanpassung
Vorhergesagte vs. beobachtete Werte
iris %>% modelr::add_predictions(mod) %>% ggplot(aes(x = pred, y = Petal.Length)) + geom_point() + geom_abline()
Weitere Diagramme zur Modellvalidierung
Residuen vs. jedes X
iris2 <- iris %>% modelr::add_residuals(mod) p <- ggplot() + geom_hline(yintercept = 0, color = "red") p1 <- p + geom_point(aes(x = Sepal.Length, y = resid), iris2) p2 <- p + geom_point(aes(x = Sepal.Width, y = resid), iris2) p3 <- p + geom_point(aes(x = Petal.Width, y = resid), iris2) gridExtra::grid.arrange(p1,p2,p3, nrow = 1)
→ Die Residuen sollten zufälligt verstreut sein. Ein Muster deutet auf eine Fehlspezifikation des Modells hin, wie hier zum Beispiel.
Visualisierung der linearen Regression
mit ggplot2
Verwende geom_smooth und setze das Methodenargument als Name der Modellierungsfunktion (hier: ‘lm’). Falls Du den Standardfehler nicht anzeigen willst, setze zusätzlich se = FALSE:
iris %>% ggplot(aes(x = Sepal.Length, y = Petal.Length)) + geom_point() + geom_smooth(method = "lm", se = FALSE)
Visualisierung der linearen Regression
manuell
iris %>% modelr::add_predictions(mod) %>% ggplot(aes(x = Sepal.Length)) + geom_point(aes(y = Petal.Length)) + geom_line(aes(y = pred))
Your turn …
Shiny App
Quiz: Verteilung & Modellspezifikation
Erkunde die Shiny App auf der folgenden Folie und beantworte diese zwei Fragen:
Shiny App (scrolle weiter runter)
Fragen?
Bei weiteren Fragen: saskia.otto(at)uni-hamburg.de
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License except for the borrowed and mentioned with proper source: statements.
Image on title and end slide: Section of an infrared satellite image showing the Larsen C ice shelf on the Antarctic Peninsula - USGS/NASA Landsat: A Crack of Light in the Polar Dark, Landsat 8 - TIRS, June 17, 2017 (under CC0 license)